% *****************************************************************************
%%                                                                            *
%%                                                                            *
% AUTHOR: Vigneshwar Ramakrishnan											  *	
% Date created: 29 August 2022												  *	
% CC-BY-SA-NC                                                                 *
%Description: this file is required for rotational_entropy.m               *
%******************************************************************************

% function [g_t0, area] = symmfft (time,vacf)

% Rotational Entropy Calculation

clear all;
% 
a = textread('c_rot.txt');
b = a';
a =b;
t = textread('timestep.txt');
 time = t;

t = time;


% 
% clear all;
% % 
% a = textread('c_rot.txt');
% t = textread('timestep.txt');
%  time = t/4;
% 
% t = time;
% a = vacf';
% a = C_Rot';
clear time;
kb = 1.3806503E-23;                 % boltzmann constant: m^2 kg s^-2 K^-1
T = 300;                           % Temperature: K
NA = 6.023E23;                     % Avogadro number
mwt_water = 0.018;                 % mol. weight of water: kg/mol (18 g/mol)
m_water = mwt_water/NA;            % mass of 1 water molecule: kg     
beta = 1/(kb*T);                    % used in w_s_HO calculation
h = 6.626068E-34;                  % m^2 kg s^-1
%  rho = 20.1242E27 % based on 4 Ang volume thickness (for protein-DNA)
rho = 32.927E27;
N_water = 128

x = a(:,1);
n2 = length (x); 

% Modify Sampling

n3 = 1;

m = round(n2/n3);

j = 1; time = 0;
for i = 1:m

 time(i) = t(j);
 x_mod1(i) = x(j);

 j = j + n3;
 
end

% Make the VAC symmetric

% Negative time

x_symm(1) = x_mod1(m);
for di = 2:m
    x_symm(di)=x_mod1(m-(di-1));
end

% positive time

for dj = 1:m-1
    x_symm(m+dj)=x_mod1(dj+1);
end
    
N = length(x_symm)

% Premultiply by the constants
x_mod = x_symm*(2/(kb*T))*N_water;

% Declare variables
hist_sum = 0;
trap_sum = 0;

% Frequency Calculations

dt = (time(2)-time(1))*1E-12; % seconds
Ws = 1/(dt);
freq = Ws*(((((-N+1)/2)):((N-1)/2))/N); % per second

% **** Non-Causal Fourier Transform Routine

count = 1;

for k = ((-N+1)/2):((N-1)/2)
    
% for k = 1:round(N/2)-1
     
    counter = 1;  
    
    for n = ((-N+1)/2):((N-1)/2)
       
        temp(counter) = x_mod(counter)*exp(-1i*2*pi*(k)*(n)/N);
        counter = counter+1;
        
    end

    % trapezoidal integration of temp
    
    trap_sum = 0;
    
    trap_sum = (temp(1)+temp(N))/2;
    for i = 2:N-1; 
        trap_sum = trap_sum + temp(i);
    end
           
    X(count) = trap_sum*dt;  
    
 
count=count+1;
end

mod = abs(X);
Fp = mod(1:N);% seconds
plot(freq(((N+1)/2):N)/3E10,Fp((N+1)/2:N)*3E10)

Fp_sav = (Fp((N+1)/2:N)*3E10)';
freq_sav = (freq((N+1)/2:N)/3E10)';
save('ft-1fs.txt','Fp_sav','-ascii')
save('freq-1fs.txt','freq_sav','-ascii')

% %fig = plot (freq,Fp,freq,Fp1,'.');
xlabel ('Frequency, (cm^-1)')
ylabel ('Rot Density of States , (cm)')
title ('Density of states for 1000 water')

% % g_t0 calculation
 
g_t0 = Fp((N+1)/2)

%area = trapz(freq((N+1)/2:N)/(3E10),Fp((N+1)/2:N)*(3E10))

D = g_t0*kb*T/(4*N_water*(1.0205+1.9181+2.9386)*1E-47)


% delta calculation

% rho = 20.62425*6.023*1e23; 


delta = (2*g_t0/(9*N_water))*(((pi*kb*T)/m_water)^0.5)*(rho^(1/3))*((6/pi)^(2/3)) %

 % f calculation

fluid = newton(0.2,delta)

% g_g calculation

for ci = 1:N
    
    g_g(ci) = g_t0/(1+(((pi*g_t0*freq(ci)))/(6*fluid*N_water))^2);
    
end

% g_s calculation

g_s = Fp - g_g; 

%%%%% Entropy calculation %%%%%% 

% Weighting function for Gas-like component

theta_x = 39.4411;
theta_y = 20.9844;
theta_z = 13.6970;

term1 = (pi^0.5)*(2.7183^1.5);
term2 = (T^3/(theta_x*theta_y*theta_z))^0.5;

s_R = log (((term1/2)*term2))*kb;
w_gs = s_R/(3*kb)

% Weighting for solid-like component

bhw = beta * h * freq;

for i = ((N+1)/2):N % 
    w_s_HO(i) = (bhw(i)/(exp (bhw(i)) - 1)) - log (1-exp (-bhw(i))); 
    % w_s_HO(i) = 1 - exp (-bhw(i));
end 

w_s_HO(((N+1)/2)) = 0;

% Integrating the weighted density of states

% solid component

for ai = ((N+1)/2):N
    
    solid(ai) = g_s(ai)*w_s_HO(ai);

end

solid((N+1)/2) = 0;

P_solid = trapz(freq(((N+1)/2):N),solid(((N+1)/2):N));

for bi = 1:N
    
    gas(bi) = g_g(bi)*w_gs;
end

P_gas = trapz(freq(((N+1)/2):N),gas(((N+1)/2):N));

Entropy = P_solid+P_gas; % For N_water water molecules

Finalvalue = Entropy*kb*NA/N_water % J mol^-1 K^-1

% ***** Energy Calculations *******

% Weight functions

for ki = 1:N
    w_E_HO(ki) = (bhw(ki)/2)+(bhw(ki)/(exp(bhw(ki))-1));
end

w_E_HO((N+1)/2) = 1;

for li = ((N+1)/2):N
    
    E_solid(li) = g_s(li)*w_E_HO(li);

end

E_gas = g_g*0.5;

E_gas1 = trapz(freq(((N+1)/2):N),E_gas(((N+1)/2):N));
E_solid1 = trapz(freq(((N+1)/2):N),E_solid(((N+1)/2):N));

Energy = E_gas1+E_solid1;

Final_energy = (Energy*NA*kb*T)/N_water;




